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We develop worldline numerical methods, which combine string-inspired with Monte-Carlo tech- 
niques, for the computation of the vacuum polarization tensor in inhomogeneous background fields 
for scalar QED. The algorithm satisfies the Ward identity exactly and operates on the level of renor- 
malized quantities. We use the algorithm to study for the first time light propagation in a spatially 
varying magnetic field. Whereas a local derivative expansion applies to the limit of small variations 
compared to the Compton wavelength, the case of a strongly varying field can be approximated 
by a derivative expansion for the averaged field. For rapidly varying fields, the vacuum-magnetic 
refractive indices can exhibit a non-monotonic dependence on the local field strength. This behavior 
can provide a natural limit on the self-focussing property of the quantum vacuum. 
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I. INTRODUCTION 

The study of vacuum polarization has been one of the 
cornerstones of the development of quantum field theory. 
It gives access to charge and field strength renormaliza- 
tion and provides for an understanding of fluctuation- 
induced modifications of Coulomb's law in terms of the 
Uehling potential which contributes to the Lamb shift. 

Vacuum polarization in strong magnetic fields has lead 
to the prediction of nonlinear optical properties of the 
quantum vacuum [Ij |2l |4H6] , for a review, see . More 
specifically, the magnetized quantum vacuum is bire- 
fringent for low-energy photons and exhibits dichroism 
above the threshold for pair production [3 ^ . Past and 
present experiments, such as BFRT [10], PVLAS [TT] . 
BMV IF, Q&A [Tg, OSQAR [H], have been working 
on the discovery of these elementary properties of the 
quantum vacuum by means of strong macroscopic mag- 
netic fields. In addition to being a fundamental test of 
QED, these experiments can provide for the strong labo- 
ratory bounds on the existence of hypothetical particles 
such as minicharged or axion-like particles [T5l - fT8] . 

As first suggested in Q^, also high-intensity lasers in 
combination with latest methods in X-ray polarimetry 
PU] may contribute to the search and discovery of these 
nonlinear properties of the quantum vacuum. Various 
laser-induced quantum- vacuum phenomena have been in- 
vestigated recently [nHHIMli, for reviews see [3DH32- 
They also provide for a powerful probe for new elemen- 
tary particles complementary to accelerators [331 - I35] or 
supplement searches for a noncommutative structure of 
spacetime [55] . 

Theoretical analyses of the vacuum polarization tensor 
so far dealt with the idealized limit of constant, homoge- 
neous fields [51[51[51 l57Hi5] which is appropriate as long as 
the scale of field variation is much larger than the Comp- 
ton wavelength in QED. Whereas this is well satisfied for 
experiments with large dipole fields searching for effects 
from fluctuations of standard model particles, this as- 
sumption can be violated for secondary fields (e.g., from 
higher harmonics) from strong lasers or for fluctuations 
involving hypothetical very light particles. A proper in- 



terpretation of future experiments will therefore require 
the knowledge of the vacuum polarization tensor in gen- 
eral inhomogeneous fields. 

Recent years have witnessed a variety of advances for 
strong-field calculations in inhomogeneous fields. So far, 
new techniques have concentrated mainly on the effec- 
tive action or effective Lagrangian in strong-fields as a 
primary quantity of interest. In addition to exact so- 
lutions [46l|47], semiclassical [48], instanton techniques, 
and quantum kinetic equations have been developed and 
applied to pair production in inhomogeneous fields, i.e., 
the imaginary part of the action, as reviewed in |31j . As 
a general-purpose numerical method, the combination of 
the worldline formalism [151 lil] with Monte-Carlo path 
integration techniques has proved successful in many in- 
stances [551156] . As one advantage, worldline numerics is 
capable of providing local information about the effects 
of fluctuations, such as energy or action densities, local 
production rates etc. Local quantities are particularly 
indicative for the non-local features of quantum-field the- 
ory. 

The present work is devoted to generalizing these 
worldline methods to the vacuum polarization tensor as 
the lowest nontrivial correlation function of QED. As the 
effective action is the generating functional for IPI cor- 
relation functions, the generalization at first sight seems 
straightforward. However, the relation is provided by 
a functional differentiation which is difficult to imple- 
ment reliably in a numerical method. In addition to 
being a powerful numerical method, the success of nu- 
merical worldline techniques also relies on the fact that 
the formulation is very close to analytical calculations. 
In fact, closed- form worldline expressions for correlation 
functions to arbitrarily high order (master formulas) can 
be derived within perturbation theory [52] . 

In this work, we demonstrate that the corresponding 
formula for the vacuum polarization tensor in an external 
field can be used to develop a worldline numerical algo- 
rithm that (i) satisfies the Ward identities at any level of 
discretization, and (ii) operates on the level of renormal- 
ized quantities so that only finite quantities are subject 
to numerical evaluation. For numerical simplicity, we 



perform all computations within scalar QED - the gen- 
eralization to spinor QED within the worldline approach 
is straightforward [511 [57] . 

As a concrete example, we concentrate on the evalua- 
tion of the refractive indices of the magnetized quantum 
vacuum deduced from the vacuum polarization tensor 
and their dependence on spatial variations. This is not 
only a crucial observable in birefringence experiments, 
but can describe a particular sensitivity of optical ob- 
servables to nonlocal features of the quantum vacuum. 
For instance, as the vacuum refractive index in the ho- 
mogeneous field approximation increases with the field 
strength, the quantum vacuum has a self-focussing prop- 
erty [5 8) : photons are drawn into regions of higher field 
strength, in turn amplifying the field strength even more. 
Our results provide for first indications that the refractive 
index can depend non-monotonically on the field strength 
in regions of large spatial variations, hence providing for 
a natural mechanism to limit this self-focussing property 
of the quantum vacuum. 

This article is organized as follows: in Sect. |Tlj we 
summarize the approach to the polarization tensor on 
the worldline. We pay special attention to spacetime in- 
homogeneities and develop a numerical algorithm which 
satisfies the Ward identity at any level of discretiza- 
tion. Sect, mil is devoted to benchmark tests in the 
form of comparisons with analytically known vacuum and 
constant-field cases. In Sect. |IV[ we present new results 
for the polarization tensor in an inhomogeneous magnetic 



field. 



II. 



VACUUM POLARIZATION TENSOR ON 
THE WORLDLINE 



Let us start with the worldline representation of the 
one-loop effective action of scalar QED in D Euclidean 
spacetime dimensions |52j . 
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dT 
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Vx e- ^0 d- i'^+^^^^) ^ 



where m denotes the (scalar) electron mass. A tran- 
sition to Minkowski-valued quantities will be discussed 
below. Here, the path integral is normalized to unity 
for vanishing field .4 = 0. Ultraviolet divergencies can 
be regularized with the help of the propertime T inte- 
gral, e.g., by replacing the lower integration limit T = 
by Tmin — l/A^ (propertime cutoff), or by dimensional 
or zeta function regularization. The corresponding effec- 
tive action for spinor QED looks very similar, addition- 
ally containing a spin-field coupling [52] . This effective 
action is the generating functional for all IPI correla- 
tion functions, which can be deduced from T by func- 
tional differentiation with respect to the gauge field A^. 
Alternatively, we can expand the gauge field in terms 
of a background field A^ and a sum over plane waves. 
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The second order is relevant 
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for the vacuum polarization tensor. 



r = i-te) 



dT 



tip 



T (47rT)^/2 

x{0)=x{T) 

r(2)M-[fc,,fc2;A]e2,-fO(e3). 



T T 

Vx e-^^'^^ T \ f dn I dT2 iieie*'=i^ii;2e2e*''=='='e-*^^'^^'^ \ +0{e') 

.0 



The plane wave basis also implements the transition to momentum space. At first sight, this appears less efficient, as 
the worldline integrals live in position space; an evaluation of the 2-point correlator in position space thus seems to 
be much more straightforward. However, it turns out that the position space representation naturally involves two 
path integrals (one for each internal propagator), whereas the momentum space formulation boils down to one path 
integral and thus is is numerically less expensive. The desired vacuum polarization tensor can be extracted as (part 
of) the coefficient of the polarization vectors e'^ , 
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The worldline integral can be decomposed into a path integral over worldlines with a common center of mass xqm 
and the spacetime integration over this center of masaj 
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In the case of a homogeneous background or in the vacuum case, there is no additional xqm dependence apart from 
the plane waves, such that the a;cM integration can immediately be performed, yielding a 5 function that implements 
momentum conservation, see below. In the general case, we have to be more careful and introduce a local vacuum 
polarization tensor tt'^'^ that depends on xcm, 

r(2)M-[fc^, fca; ^] - J d^xcM e^('^^+'=^)^™7r''''[fci, fc2, xcm; A], (1) 

where the local polarization tensor is 
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(2) 



In the homogeneous or vacuum case, Eq. ([I]) boils down 
to 

rf^V-Ifci^fca;^] = (27r)^5(^)(fci + fc2)n'^''[fc;A], (3) 

where Il^'^[k; A] is the standard vacuum polarization ten- 
sor in a homogeneous field, and k = fci = —k2- Mo- 
mentum conservation is automatically implemented, and 
local and global descriptions are identical. 

U(l) gauge symmetry imposes a constraint on the po- 
larization tensor in the form of the Ward identity. In its 
local form, the polarization tensor has to satisfy 

ki,tj.'^'"'[ki,k2,xcM;A] ^ TT''''[ki,k2,xcM;A]k2.^ =0. 

(4) 
In the worldline representation ([T]), the Ward identity be- 
comes obvious, as the r^ integrands turn into total deriva- 
tive upon contraction with fc^. 
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where the last equality holds as the worldlines are closed, 
Xi(0) = Xi{T). However, in any discretized numerical ap- 
proach, this identity is difficult to realize as there is no 
Leibniz rule for latticized derivatives. Nevertheless, the 



Ward identity can be exactly maintained in the numeri- 
cal worldline algorithm due to the following observation. 
Consider the parameter integral combination /^i^^ q(,_ 
curring in Eq. M, 
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The mixed terms arising from an expansion of the last 
equation vanish in the continuum by virtue of Eq. ([5|. 
The same holds for the product of the second terms in 
parentheses. This representation explicitly shows that 
the tensor structure of /'^^^^ is identical to its contraction 
with two corresponding transversal projection operators 

H 

a property which also holds for the local polarization ten- 



^ Different prescriptions for such a decomposition can equally well 
be used 52, 57). 



(8) 



by virtue of which the polarization tensor satisfies the 
Ward identity manifestly. In other words, using the iden- 
tity ([7]), the local polarization tensor yields a form 
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which satisfies the Ward identity manifestly also upon 
discretization of the worldline on a propertime lattice. 
This discretization can now proceed in the standard way 
[551 . First, we rescale the worldlines. 



(10) 



x(T) = Vry(i), ,r^Tt, 
such that the kinetic term 
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dty{tY , (11) 



serving as the probability distribution of the Monte Carlo 
configurations, becomes independent of the propertime 
T. Each worldline y{t) is then represented by a set of 
N points per loop (ppl), yi — y{ti), where ti = i/N and 
i = l,...,Af. The worldline path integral then turns 
into an expectation value with respect to an ensemble of 
a^CM centered worldlines {yi}., where £ = l,...,nL and 
riL denotes the number of worldlines in the ensemble. 
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The expectation value is normalized to (1) — 1. There 
are powerful algorithms available to generate the ensem- 
ble ab initio [SB] [52] • In the present work, we use the 
V loop algorithm [5B]. Together with the rescaling (10 1, 
the local polarization tensor then reads 



7r^''[fci,/c2,a;cM;^] 

00 

{-ief r dT 
- (4^)D/2 J Tfne 




(13) 
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This representation of the local polarization tensor serves 
as the master formula for the construction of our algo- 
rithm and its application in the following sections. 



III. BENCHMARK TESTS 

A. Vacuum polarization tensor at zero field 

At Afj_ = 0, any dependence on the spacetime coordi- 
nate drops out. Homogeneity ensures that the local po- 
larization tensor obeys 4-moinentum conservation, such 
that Eq. ((sl) applies. The worldline expression of the un- 
renormalized polarization tensor then reads 



Ut'^ik] 



(-ie)2 r dT 
(47r)^/2 J WJ^ 



e-"" ' PJ^''{k)P^\k) (14) 
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For a benchmark test in £> = 4 spacetime dimensions, we 
need to renormalize Eq. ( |14[ ) and then perform a com- 
parison with the analytically well-known result for scalar 
QED [H] 



W^ik) = 



(47r)2 



'-PJ^-'ik) 
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(15) 
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where the last term in the curly brackets corresponds to 
the counter-term from charge renormalization. 

Renormalization is not only an important conceptual 
issue, but also needs to be taken care of as a matter of 
practice, as it is advisable to perform the numerics only 
for finite renormalized quantities. As the UV divergen- 
cies associated with Eq. ( 14 1 occur for small propertimes 
T, the divergencies can be analyzed and taken care of 
on the level of the propertime integrand which is always 
finite. Only after the counter terms are subtracted, we 
perform the propertime integration which then is per- 
fectly finite. Even though propertime regularization is 
most natural for our formalism, our result can be con- 
nected to any other scheme, as the counter terms are 
known analytically from the small propertime expansion 
of the worldline integral. 



Still the subtraction is not trivial, as subtracting the 
analytically known counter term (1/3) from the numer- 
ically evaluated propertime integrand in Eq. ( 14 ) would 



lead to an indefinite result with infinitely large error bars. 
This is because the small propertime behavior of the nu- 
merical worldline expression is equal to (1/3) only within 
numerical precision. 

The solution to this problem has been provided by [S3] 
for effective action computations: the existence of the 
heat-kernel (small propertime) expansion of the world- 
line expectation value at least in an asymptotic sense 
allows us to fit the numerically obtained worldline ex- 
pression to a polynomial in T. Subtracting the constant 
piece (= (1/3) within numerical errors) then corresponds 
to charge renormalization. For improving the stability 
of the numerical result, it is advisable to fit to a higher- 
order polynomial (the coefficients of which can also easily 
be worked out analytically with the heat-kernel expan- 
sion). The propertime integrand of tt^'' is then evaluated 
with the pure numerical result for values of T larger than 
a scale Tdg and with the fit for T < Ty^q . The scale Tdg 
is dynamically generated from the condition that both 
expressions for the propertime integrand should have er- 
ror bars with the same size at Tdg- 

In Fig. [l] the numerical and analytical results of the 
propertime integrand evaluated in Euclidean space for a 
momentum vector k^i = (1,1,1,1) and Lorentz indices 
chosen in the 11-direction is shown. Here, the units are 
chosen such that each momentum component /c^ in ar- 
bitrary inverse length units L~^ has the value k^ — 1 
for /i = 1,2,3,4. In the upper panel, we have set the 
electron mass m = Q whereas the lower panel depicts 
the cases m = 1 and m = 2 (using the same arbitrary 
length units). The latter show a characteristic exponen- 
tial drop-off for large propertimes arising from the e~™ '^ 
factor. In all cases, the numerical results represent a very 
satisfactory approximation to the analytical results |42j . 

In Fig. [2J the diagonal (upper panel) and off-diagonal 
(lower panel) components of the vacuum polarization ten- 
sor in Euclidean space are shown as a function of the 
mass parameter m again for the case k^ = (1,1,1,1). 
The good agreement between analytical and numerical 
results for a wide range of mass values demonstrates that 
our method is capable of computing perturbative corre- 
lation functions with worldline numerics. 

As described above, we have performed the renormal- 
ization by polynomially fitting the propertime integrand 
and successively subtracting the counterterm correspond- 
ing to charge renormalization. More precisely, the term 
in curly brackets in Eq. ( 15 ) is fitted, for instance, to 
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Figure 1. Comparison of the analytical and numerical re 
suits of the propertime integrands of Eqs 
(1 



(14} and ([15|) for 
(in arbitrary inverse length units L~ 



kfj, = (1,1,1,1) (in arbitrary inverse length units L~ and 
Lorentz indices chosen to point into the 11-direction for dif- 
ferent values of the electron mass: m = (upper panel) and 
m = 1 and m = 2 (lower panel). In all cases, we have used the 
same ensemble of random worldlines. A test of the algorithm 
shows that an ensemble with A'^ppi — 1000, riL = 40000 gives 
us acceptable results with respect to both the calculation time 
and the numerical errors. 



heat-kernel expansion can as well be employed in the case 
of nonvanishing electromagnetic fields. 



IV. POLARIZATION TENSOR IN A 
HOMOGENEOUS MAGNETIC FIELD 

As is obvious in the worldline approach, the generaliza- 
tion to nonvanishing background fields is straightforward, 
by inserting the Wegner- Wilson loop 



ef^dr x^A'^{x{t))^ 



into the worldline average in Eq. (14), cf. Eq. (13), with 
x{t) = xcM + VTy{t), t — t/T. As another test of our 
method, we compute the vacuum refractive indices aris- 
ing from fluctuations in a homogeneous magnetic field. 
Due to homogeneity, the 4-inomentum of the photon is 
conserved, implying ki + k2 = 0, such that we can di- 
rectly study the polarization tensor. Writing the gauge 
potential in the form 



P = bT + cT^ 



dT'' 



(16) 



or higher-order polynomials. The algorithm can be sta- 
bilized by inserting the analytically known coefficients 
from the heat-kernel expansion, b = — fc^/30, c = A:''/420, 
d = — fc^/7560. Similar techniques and knowledge of the 



^M ^ Ip'^'i^ = -BA'^'i^ 



with the magnetic field strength B and a dimensionless 
4-diniensional tensor A'^". Choosing the magnetic field 
to point into the ei direction. A'''' in Euclidean as well 
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2 - Point Function with k = (1 , 1 , 1 , 1 )""" 




Figure 2. Comparison of the analytical and numerical results for the full polarization tensor If" in scalar QED as a function 
of the electron mass m. The upper panel depicts a generic result for the diagonal elements (/i = i/). The lower panel shows the 
same for off-diagonal elements. 



as in Minkowski space reads j^ 
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The refractive indices are the inverse of the phase veloc- 
ities of photons propagating in a magnetized quantum 
vacuum. As the magnetic field distinguishes a direction 
in space, the magnetized quantum vacuum is birefringent 
like a uniaxial crystal, featuring two polarization depen- 
dent phase velocities. 



II/-L 



= 1- 



Hi 



,_L 



k2 



= (l-At'||/±r 



(17) 



where nn j^ are the nontrivial eigenvalues of the eigen- 
modes e^ of the polarization tensor 0, H'^'^e^.y/.L = 
n||/_Le;j,||/_L satisfying the Minkowski-space photon dis- 
persion relation 



■n,. 



= 0. 



(18) 



^ In the case of a Minkowskian electric field, imaginary Euclidean 
components would have to be inserted into the worldline inte- 
grals, see f55| . 



Equations (171 and ( |18| are actually identical as the 
phase velocity is defined as v — tj/|k|. Here, we use 
the Minkowski metric g = (— , +, +, +) and parameterize 
the Minkowskian momentum as (fc^)M — {^i\M)- For 
scalar QED, the velocity shifts, ^vyj_ = 1 — i;||_l in the 
weak-field limit yield [T7] 
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(19) 



Here 9 denotes the angle spanned by the magnetic field 
B and the propagation direction k. For weak fields, the 
II mode is polarized in the plane spanned by B and k, 
whereas the _L mode is polarized orthogonal to this plane; 
more explicitly, ey ^ Afc and e± ^ Ak, with A denoting 
the dual field strength matrix. 

Even though the worldline numerical formalism is 
set up in Euclidean space, as the Monte Carlo proce- 
dure requires a positive action, extracting these light 
propagation properties obviously requires a transition 
to Minkowski space. In particular, we have to insert 
a Minkowski-valued 4-momentum vector into the world- 
line average to have access to real light-cone properties. 
We do so by choosing the Euclidean 4-momentum vec- 
tor as fc^i = (iWjk). As illustrated in the Appendix, the 
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Figure 3. Benchmark test: velocity shift prefactor a|| in a 
homogeneous magnetic field for the || mode at orthogonal in- 
cident 6 = 7r/2 and for uj = O.lm. For weak fields, ay ap- 
proaches the analytical result 1/90 (« 0.01111), see Eq. (19 1. 
Also for larger fields, ay follows the analytically known non- 
perturbative result |17| . 



algorithm remains stable at moderate frequencies, even 
though larger fluctuations require typically two orders of 
magnitude more statistics than typical Euclidean com- 
putations. 

In Fig. [3] we compare our worldline Monte-Carlo re- 
sults with those of the analytically known velocity shifts 
[T7] over a wider range of magnetic field strength (also 
exceeding the simple weak field limit). This benchmark 
test is performed for orthogonal incident 9 = 7r/2 and 
for the II mode at a frequency w = O.lm. The numerical 
results approach the velocity shift in the weak-field limit 
with very good accuracy and also yield reliable results 
for larger field strengths. 



V. POLARIZATION TENSOR IN A SPATIALLY 
INHOMOGENEOUS MAGNETIC FIELD 

Let us now explore new vacuum polarization effects 
in inhomogeneous fields, revealing the nonlocal nature of 
fluctuation-induced processes. For this, we use a mag- 
netic background pointing into, say, ei direction, con- 
sisting of a constant magnetic field B superimposed with 
a sinusoidal oscillation varying in 63 direction with am- 
plitude Bi and wavelength A^, 



B(a:3) 



, 27r 
B + Bi cos 1 T— a:^3 
As 



ei. 



(20) 



A similar electric field has already been used to ana- 
lyze the role of spatial inhomogeneities in Schwinger pair 
production [55]. This field configuration can be viewed 
as a rough approximation to a realistic strong and broad 
laser pulse in standing wave mode superposed with higher 
harmonics. 

For the worldline simulation, we use the corresponding 
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Figure 4. Sketch of the geometry of our configuration with 
an inhomogeneous magnetic field B(a;3) (green arrows and 
corrugated surface) with a constant magnetic component B 
in ei direction and a spatial variation along the 63 direction 
with field amplitude B\ also pointing into the ei direction. 
The photon (blue wiggly line) propagates perpendicular to 
the field in 62 direction. The fluctuation-induced interaction 
between the photon and the magnetic field is represented by 
a Feynman diagram (black ellipse). 



gauge potential 
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^ 27r 



27r 
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which is numerically convenient, as it depends only on 
one spatial coordinate. 

As a relevant observable, we compute the local velocity 
shift Aw (2:3) for a photon propagating along the 82 direc- 
tion from the local polarization tensor. The geometry of 
our configuration is sketched in Fig. |4] As the magnetic 
field is homogeneous in 62 direction, the photon momen- 
tum is conserved /ci + ^2 = for ki ^ k2 ^ 62. Hence, 
setting k — ki = — ^2 and inserting (21 1 into the local po- 
larization tensor ( 13 1, we can determine T:'^'^[k, xs^cm! A], 



which is diagonalized by the same polarization eigen- 
modes e^i.||/± as the in constant-field case. The local 
phase velocity shifts then are computed analogous to 
Eq. (fTTl), 



.(2^3^ 



_(a;3. 



CM, 



11,^. „,CMj , . , ,,2 

= 1 r^ = (1 - ^«||/_l(2^3,cm)) , 

(22) 
where ttii _L(a;3_cM) £^re the local eigenvalues of the polar- 
ization tensor. 

In our cornputation, we use eB = 0.2m^ and Bi = 
0.5B. Figure [5] shows a contour plot of the local velocity 
shift as a function of the oscillation wavelength Xb and 
the normalized coordinate a:3_cM/As for the ||-mode; the 
latter corresponds to the phase inside the oscillation pe- 
riod, (j) = 27rx3, cm/As. For large Xb when the field be- 
comes slowly varying with respect to the Compton wave- 
length Xbtti ^ 1, the local velocity shifts approach the 




Figure 5. Contour plot of the phase velocity shift Av for 
an incoming ||-photon at orthogonal incident with respect to 
the direction of the external field (6' = ^). The inhomoge- 
neous field has a sinusoidal variation along the 63 direction, 
cf. Eq. Eo\, with eB = 0.2m^ and Bi = 0.5B. The phase ve- 



locity is plotted horizontally versus the oscillation wavelength 
Xb and vertically versus the normalized position 0:3, cm/As of 
the photon relative to the oscillation phase. Lines of con- 
stant absolute position a;3_cM = 0.5, 1.0,2.0 in units of mass 
m for varying wavelength As are indicated as dotted, solid 
and dashed lines. 
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Figure 6. Contour plot with the same physical parameters 
as in Figure |5] but with a higher resolution in the region of 
small Xb- The y axis is now given in terms of the phase in 
units of TT, (/!>/7r — 2xcm/Xb- The straight lines are cuts at 
A = 0.15, 0.4, 0.6, which we discuss in Figure |8] 



constant field values as expected. A "locally-constant- 
field" approximation becomes reliable in this limit. By 
contrast, if the two characteristic length scales become 
similar A^m ~ 1 the oscillating structure of the mag- 
netic field starts to become washed out in the local ve- 
locity shift. The propagating photon undergoing virtual 
electron-positron loops with an inherent length scale of 
the Compton wavelength "sees" a field averaged over this 
length scale. For a rapidly oscillating field, A^m <C 1, the 
field oscillations are completely washed out and become 
invisible in the velocity shift. In our present example, 
the limiting velocity shift in this rapid-oscillation limit 
corresponds precisely to that induced by the background 
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Figure 7. Phase velocity shift A.v\\ for different positions in 
the phase of the variation </; = 2'KXz,cyilXB as a function of 
the variation length As. The curves correspond to horizon- 
tal slices of the contour plot Figs. [SJGl at the ordinate values 
(ji/vr = 22;cm/As = 0,0.5, 1. For large As, the velocity shifts 
approach their constant-field limits with a clear ordering from 
large to small background field from top to bottom (blue dots 
at field maximum to red stars at minimum field). By contrast, 
this ordering is modified at small As in various patterns de- 
pending on the value of As (see inlay). The straight lines in 
the inlay at As — 0.25, 0.65 mark the turning points. 



field i?. Aw ^ B . This is in line with the interpreta- 
tion that the averaging arises from the nonlocal nature 
of the fiuctuations. For instance, naively averaging over 
the constant-field velocity shift, depending quadratically 
on i?, would give a different (and wrong) averaging result, 
Ai>7^(# + iB2). 

An interesting parameter regime occurs at small Xb-, 
i.e., for rapidly varying background field, see Fig. [61 In 
Fig. [7) we show horizontal slices of the contour plot Fig. [6] 
at the ordinate values Xs^cm/Ab — Oj 0.25, 0.5, corre- 
sponding to different positions in the phase of the vari- 
ation (/)/7r = 2a::3_cM/A_B = 0,0.5, 1, as a function of the 
variation length A^. For large A^, the velocity shifts ap- 
proach their constant-field limits with a clear ordering 
from large to small background field from top to bot- 
tom. This is in accordance with expectations from a 
locally-constant-field approximation becoming applicable 
for large \b- 

By contrast, this ordering is lost at small As, where the 
curves show a characteristic oscillation pattern, see inlay 
of Fig. [7] Depending on the value of A^ , the velocity shift 
in a local minimum of the field (red stars, (/)/7r = 1) can 
become larger than that in a local maximum (blue dots, 
(/)/7r = 0), cf. inlay of Fig. [t] at around Xb — 0.65. We 
interpret this phenomenon as a consequence of the local 
averaging property of the quantum fluctuations on scales 



of the Compton wavelength m. In this way, the local 
velocity shift in a minimum of the field can receive dom- 
inant contributions from the nearby maxima if they are 
significantly probed by the quantum fluctuations on the 
scale 1/m. Conversely, the local velocity shift in a maxi- 
mum of the field can receive dominant contributions from 
the nearby minima. This can lead to an inversion of the 
hierarchy of the velocity shifts with respect to the local 
background field. Our data is compatible with further os- 
cillations setting in at even smaller values of Xb — 0.25, 
which would correspond to further minima or maxima 
entering the local fluctuation average. 

This becomes visible in vertical slices of the contour 
plot Fig. |6] at the abscissa values Xb = 0.15,0.4,0.6, 
and 2 as a function of the phase of the field variation 
(j) = 2TTX3fiM/XB- In Fig. pi we show the varying part of 
the velocity shift normalized by the maximum oscillation 
amplitude A. We observe the expected 27r-periodicity 
of the phase velocity with respect to (p. Most impor- 
tantly, the oscillation is shifted by an offset (j>o ~ tt for 
small wavelengths, e.g., Xb = 0.4, 0.6, compared to larger 
As = 2 where the velocity shift tends to approach the 
locally-constant-field limit. This implies that velocity 
shift minima occur at field maxima and vice versa in 
the range 0.25 ^ As < 0.65. For even more rapid field 
oscillations A^ < 0.25, our data is compatible with the 
velocity shift being in phase with the external field again. 
This goes hand in hand with our interpretation in terms 
of fiuctuation averages. The full Monte Carlo data can 
be parametrized by a simple fit. 



f{x) = Acosi 27T 



Ar 



A Wo. 



(23) 



The fitting results are given in Tab. IT] 

A similar phenomenon had already been observed for 
the case of pair production in inhomogeneous fields |55| ; 
however, due to the exponential dependence of pair pro- 
duction on the background field, this averaging phe- 
nomenon was much more pronounced in this case. In 
fact, the amplitude of the oscillations is on the order of 
our statistical error bars for the velocity shift data. How- 
ever, fitting the data at fixed Xb to sinusoidal fit func- 
tions translates into a correspondingly large error for the 
amplitude, but a significantly small error for the phase of 
the oscillation, see Tabjlj The size of our statistical errors 
can also be estimated from the deviations from symme- 
try about the horizontal x^^qm/^b — 0.5 axis in Fig. [5] 
or from anti-symmetry about the horizontal </> = 0.5 axis 
in Fig. [6J This latter anti-symmtry also guarantees that 
the exact result for the phase velocity shift at = 0.5 
should not depend on A^; the slight dependence of the 
Monte Carlo data for the = 0.5 curve (green diamonds) 
on As, hence is a measure for our statistical error rl 
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Figure 8. Normalized phase velocity shift Auy for different 
field inhomogeneities As = 0.15,0.40,0.6,2.0 as a function 
of the phase of the field variation (J)/tt — 2x3^cm/^b- The 
curves correspond to vertical slices of the contour plot Fig. |6] 
at the given abscissa values for A^. The picture shows the 
expected 27r-periodicity of the phase velocity with respect to 
(j). For large As = 2 (light blue stars), the velocity shift is in 
phase with the background field ~ cos (f> as expected from a 
locally-constant field approximation. For As = 0.4, 0.6 (green 
diamongs, red triangles), the velocity shift oscillates out of 
phase ~ cos{(j) — cjto with (j)o ~ vr. For even smaller variation 
length As = 0.15, the oscillation is approximately in phase 
again ~ cos <j). 



^ Of course, these (anti-)symmetries could be implemented explic- 
itly in the Monte Carlo computation by generating correspond- 
ingly symmetric ensembles. 



As these velocity shifts correspond to shifts of the local 
refractive indices this new phenomenon has a direct con- 
sequence for the self-focussing property of the quantum 
vacuunr [58]: in the locally-constant-field limit (a pure 
Heisenberg-Euler-type calculation), the refractive index 
increases with increasing field strength. This implies 
that photons are dragged into local maxima of the field 
strength. This even enhances the field strength at local 
maxima, thus giving rise to self-focussing properties. Our 
observation in turn predicts that this self-focussing natu- 
rally terminates on the scale of the Compton wavelength. 
If field maxima are self-focussed down to a critical scale 
Acr.i, field maxima with nearby minima can become lo- 
cal minima of the velocity shift (and thus minima of the 
local refractive index) such that the quantum vacuum 
becomes defocussing again. From the inlay of Fig. [T] we 
estimate this critical scale to be near Acr,i — 0.65 in units 
of the Compton wavelength. This critical scale provides 
a natural limit to the self-focussing property of the quan- 
tum vacuum. As our data is compatible with an in-phase 
dependence of the refractive index on the field inhomo- 
geneities for Xb < 0.25, the quantum vacuum may be- 
come self-focussing again on this shorter variation scale. 
But this regime is also expected to terminate at another 
critical scale Acr,2 where the velocity shift may run out of 
phase again. Our data is compatible with Acr,2 ^ 0.05. 
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variation length \b 


amplitude A 


phase shift (/io/ti" 


constant velocity shift Auq 


0.15 


7.36-10-^±2.5-10-^ 


0.31±0.30 


2.62-10-''±3.4-10-* 


0.40 


1.61-10-*±4.8-10-® 


1.09±5.4-10~^ 


2.62-10-''±5.1-10-* 


0.60 


8.66-10"''±3.4-10~* 


1.13±6.8-10-' 


2.62-10-^±5.M0"* 


2.00 


8.52-10"^±2.1-10"^ 


0.02±3.0-10"^ 


2.74-10"®±3.8-10"* 



Table I. Fit parameters deduced from the numerical results of Fig. [s] using a least-squares fit to Eq. (23 1. Starting from large 
A_B (slowly varying field), the velocity shift is in phase with the external field, (f)o — 0. For As = 0.4, 0.6, we observe a jump of 
the phase by ^o/tt — 1. For an even smaller variation scale As — 0.15 the phase 4>o goes back to values compatible with zero 
(the deviations from zero arise from the comparatively large signal to noise ratio of our data for small As ; this is also reflected 
in the large error bars for the amplitude A for smaller As). 



VI. CONCLUSIONS 

Based on the successful virorldline approach to per- 
turbative correlation functions, we have developed nu- 
merical Monte-Carlo techniques for the computation of 
the vacuum polarization tensor in inhomogeneous back- 
ground fields for scalar QED. These techniques generalize 
earlier methods which have been frequently applied to ef- 
fective action or quantum energy computations. The new 
challenge in the case of correlation functions is the ap- 
pearance of further scales provided by the incoming and 
outgoing momenta. 

The stability of the numerical algorithm also originates 
in the fact that it satisfies the Ward identity exactly 
and operates on the level of renormalized quantities. We 
have explicitly demonstrated that the algorithm can also 
be used to determine correlation functions as a function 
of Minkowski-valued momenta and fields, even though 
stability is expected to become an issue for increasing 
Minkowski momenta or dominating electric field compo- 
nents. 

We have verified our algorithm with the analytically 
known cases of the vacuum polarization tensor for off- 
shell momenta and the polarization tensor in homoge- 
neous fields using the magnetically induced light-cone 
deformations as an observable. In these cases, the al- 
gorithm is capable to reach a precision on the percent 
level at moderate numerical cost. 

Furthermore, we have studied for the first time light 
propagation in a spatially varying magnetic field. For 
small variations of the field compared to the Compton 
wavelength, the local derivative expansion (or locally- 
constant-field approximation) is well applicable as ex- 
pected, such that the vacuum polarization tensor quickly 
approaches the constant-field limit. 

For rapidly varying fields, the vacuum-magnetic refrac- 
tive indices can exhibit a non-monotonic dependence on 
the local field strength. This new behavior can geomet- 
rically be understood in the worldline picture, as the 
worldlines and their spatial extent probes the nonlocal 
structure of quantum field theory. Local values of the re- 
fractive indices can receive dominant contributions from 
nearby maxima or minima of the field strength. This 
inherent averaging mechanism induces a smearing and 
even non-monotonical features of the refractive indices. 



For the properties of light propagation, this can provide 
a natural limit on the self-focussing property of the quan- 
tum vacuum. 

Our present study represents a first step into the 
largely unknown territory of quantum correlation func- 
tions in inhomogeneous fields. Even though we have 
concentrated on the two-point function in the present 
work, we expect that our algorithmic strategy can rather 
straightforwardly be generalized to higher-order correla- 
tion functions. Also a generalization to spinor QED is in 
principle straightforward and merely requires the inclu- 
sion of the worldline spin factor involving a numerically 
moderately expensive path-ordering prescription. 
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Appendix A: Validity control of the numerical 
algorithm 

The worldline Monte Carlo method has proven its 
efficiency and accuracy in many examples in the con- 
text of effective action and quantum energy computa- 
tions. Generically, the convergence is very satisfactory 
and scales with a typical Monte Carlo 1/y/nZ depen- 
dence, where riL is the number of configurations, i.e., 
worldlines in this case. Precision with an error on the 1% 
level can be achieved with moderate numerical cost. In 
the present case, it is worthwhile to critically re-examine 
the quality and efficiency of the algorithm, as the cal- 
culation of correlation functions goes along with further 
technical requirements. Most prominently, the physical 
observables need to be computed with Minkowski-valued 
momenta which technically is a potential source of nu- 
merical instability. This is because the Euclidean phase 
factors '^ e'x.'p{ikf^Xfj) receive real exponential contribu- 
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♦ constant field (minkowskian vector) 

• constant field (euclidean vector) 

* inhomogeneous field witfi ?^ = 1 .0 

* inhomogeneous field with X = 7.0 
function ~ 1 / n"^ 



Number of Loops n 

Figure 9. Relative error (jackknife estimate) for the veloc- 
ity shift Av in percent as a function of wl for a Euclidean 
vs. Minkowskian calculation including its dependence on the 
field inhomogeneity. We have used a worldline ensemble with 
A^ = 1000 points per loop. The plot shows a comparison be- 
tween the Euclidean constant-field case with eB/rn^ = 0.2 
and a Minkowskian calculation for the same constant-field 
case (blue diamonds). Two Minkowskian calculations for the 
inhomogeneous field of Sect. [V] with wavelengths As — 1.0 
(red stars) and As — 7.0 (light blue triangles) are also shown. 
Because of the stronger fluctuations in the Minkowskian case, 
all Minkowskian calculations show larger relative errors. The 
dependence of the error on the background inhomogeneity is 
comparatively minor. For illustration, also a line ~ l/y^nE is 
shown reflecting the expected error depletion for Monte Carlo 
calculations. 



tions for Minkowskian momenta ko — ^ ico. Also, the 
frequency and spatial momentum dependence introduces 
further scales which can interfere with the scale of inho- 



mogeneity. 

All error estimates in this work are based on the Jack- 
knife method. We have checked explicitly, that this error 
estimate using the same random number seed yields re- 
sults equivalent to an error estimate derived from a set of 
ensembles created with different random number seeds. 

In Fig. |9l we compare the relative error for the veloc- 
ity shift Awii in percent as a function of the number of 
worldline configurations tt-l for various parameters. The 
smallest error is observed for a purely Euclidean constant 
field calculation (green dots) . This type of calculation is 
closest to conventional effective action and quantum en- 
ergy computations; however, in the present case it does 
not describe a physical observable. The corresponding 
physical Minkowskian calculation for the constant-field 
case (blue diamonds) shows an error increase of roughly 
an order of magnitude. This implies that an error on 
the level of a typical Euclidean calculation requires two 
orders of magnitude more statistics. 

The interference of the error with the external scales is 
also illustrated in Fig. |9J where we determine the relative 
error for the velocity shift Av in percent as a function of 
TT-L and its dependence on the field inhomogeneity for two 
different values Xb parameterizing the field inhomogene- 
ity. Even though we observe some dependence on the 
field inhomogeneity leading to slightly larger errors, the 
main effect on the error clearly arises from the necessity 
to perform Minkowski-valued calculations. 

Fig.lolalso depicts a straight line exhibiting a ^ 1/^nL 
dependence in this double-log plot. This indicates that 
the numerical error decreases with increasing number of 
worldhnes riL as ~ l/y/uh, as expected. Even for these 
perturbatively small Minkowski-valued quantities, errors 
below the 10% level are achievable at managable numer- 
ical cost. 
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